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Abstract. A key issue in dimension reduction of dissipative dynamical systems with spectral 
gaps is the identification of slow invariant manifolds. We present theoretical and numerical results 
for a variational approach to the problem of computing such manifolds for kinetic models using 
trajectory optimization. The corresponding objective functional reflects a variational principle that 
characterizes trajectories on, respectively near, slow invariant manifolds. For a two-dimensional linear 
system and a common nonlinear test problem we show analytically that the variational aijproach 
asymptotically identifies the exact slow invariant manifold in the limit of both an infinite time horizon 
of the variational problem with fixed spectral gap and infinite spectral gap with a fixed finite time 
horizon. Numerical results for the linear and nonlinear model problems as well as a more realistic 
higher-dimensional chemical reaction mechanism are presented. 
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1. Introduction. In dissipative ordinary differential equation systems modeling 
chemical reaction kinetics the phase flow generally causes anisotropic volume contrac- 
tion due to multiple time scales with spectral gaps. This leads to a bundling of 
trajectories near "invariant manifolds of slow motion" of successively lower dimension 
during time evolution. Model reduction methods exploit this for simplifying the un- 
derlying ordinary differential equation models via time scale separation into fast and 
slow modes and eliminating the fast modes by enslaving them to the slow ones as a 
graph of a function which deflnes the slow invariant (attracting) manifold (SIM). 

Early model reduction approaches in chemical kinetics like the quasi steady-state 
and partial equilibrium assumption [32] have been performed "by hand" , modern nu- 
merical approaches are supposed to automatically compute a reduced model without 
need for detailed expert knowledge of chemical kinetics by the user. Many of these 
techniques are based on an explicit time-scale analysis of the underlying ordinary 
differential equation (ODE) system. 

Among those methods that became popular in applications are the intrinsic low 
dimensional manifold (ILDM) method [21] and recent extensions of its main ideas, e.g. 
the global quasi-hnearization (GQL) [4], computational singular perturbation (CSP) 
[15, 16], Eraser's algorithm [6, 9, 23], the method of invariant grids [5, 11, 12], the 
constrained runs algorithm [10, 33], rate-controlled constrained equilibrium (RCCE) 
[14], the invariant constrained equilibrium edge preimage curve (ICE-PIC) method 
[27, 28], flamelet-generated manifolds [7, 30], and finite time Lyapunov exponents 
[22]. For a comprehensive overview see e.g. [11] and references therein. 

Reaction trajectories in phase space that are solutions of an ODE system x{t) = 
f{x{t)),x{0) = xo,f € C°°, describing chemical kinetics are uniquely determined 
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by their initial values and the corresponding orbits bear global information about 

phase space structure. Based on Lebicdz' idea to search for an extremum principle 
that distinguishes trajectories on or near slow attracting manifolds, an optimization 
approach for computing such trajectories has been applied in [17, 18, 26]. In [19] the 
authors propose and discuss various geometrically motivated optimization criteria for 
the formulation of a suitable extremum principle and present numerical results for 
several applications. 

The present work systematically analyzes a variational formulation of the prob- 
lem to compute slow invariant manifolds and its potential for identifying the correct 
manifold for linear and nonlinear test problems in two-dimensional phase space. We 
analytically prove the c;orrect identification of the slow eigenspace and SIM respec- 
tively in the limit of infinite-time horizon of the variational problem and derive an 
error quantification as a function of spectral gap and finite time horizon length. In 
addition, we provide corresponding numerical results confirming the theoretical pre- 
diction. 

2. Variational Problem. We consider autonomous ODE systems of the form 
X = f{x) modeling chemical reaction kinetics that have a stable fixed point corre- 
sponding to chemical equilibrium. The basic idea of our approach is the formulation 
of a variational principle that captures essential properties of a slow invariant man- 
ifold (SIM). We propose an appropriate characterization of maximum "slowness" in 
terms of an integral over suitably defined curvature (velocity change) of trajectories 
measured in the Euclidean norm. The SIM is generally characterized by the property 
that all trajectories in its neighborhood converge faster to the manifold than to the 
attractor, the chemical equilibrium point. Adrover et al. [1] recently argued that this 
might be interpreted as a ratio r > 1 of the local stretching (contraction) rate of 
vectors orthogonal to the manifold compared to those tangent to the manifold. This 
point of view comes close to our reasoning on the basis of a variational principle. 

2.1. Trajectory-Beised Optimization Approach. The variational problem 
can be formulated as 

min / ^(x(t)) dt (2.1a) 

subject to 



dx{t) 
dt 



f{x{t)) (2.1b) 
= g{x{t,)) (2.1c) 

Xj{t^) = x]', je /fixed, (2.1d) 



with to ^ ^ tf. The variable x = {xi)"^i denotes the state vector and /fixed is 
an index set that contains the indices of state variables (denoted as reaction progress 
variables in chemical kinetics) with fixed values at fixed time chosen to parame- 
terize the reduced model, i.e. the slow attracting manifold to be computed. Thus, 
those state variables representing the actual degrees of freedom within the optimiza- 
tion problem are Xj{t^.), j ^ /fixed- The process of determining .x'',.?' ^ /fixed from 
a;** , j e /fixed is known as species reconstruction in chemical kinetics and represents 
a function mapping the reaction progress variables to the full species composition 
by determining a point on the slow attracting manifold. The system dynamics (e.g. 
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chemical kinetics determined by the reaction mechanism) are described by (2.1b) and 
enter the optimization problem as equality constraints. Hence an optimal solution 
of (2.1) always satisfies the system dynamics of the full ODE system and therefore 
represents a solution trajectory of (2.1b). Additional constraints (e.g. chemical el- 
ement mass conservation relations in the case of chemical kinetics that have to be 
obeyed due to the law of mass conservation) are collected in the function g in (2.1c). 
The state variables chosen as parameterization of the reduced model (slow invariant 
manifold) are fixed via the equality constraint (2.1(1) at t*. The objective functional 
<^{x{t)) in (2.1a) characterizes the variational principle that will be discussed in the 
next section. 

2.2. Optimization Criterion. In [19, 26] 



is proposed as a suitable criterion with Jf{x) being the Jacobian of the right hand 
side / evaluated at x{t) and || • ||2 denoting the Euclidean norm. 

The term Jf{x) f represents the rate of change of reaction velocity in its own 
direction along a trajectory and can be interpreted as a specific definition of curvature 
in time parameterization of the curve 



The minimization of the time integral over <I> in (2.1a) incorporates the "maximum 
slowness" issue in terms of an average over suitably measured local curvature of a 
trajectory. 

For further analytical and numerical investigation of the variational formulation 
we consider a slight modification of the objective functional 



2.3. Forward and Reverse Mode. In previous publications [17, 18, 19] the 
general optimization problem (2.1) is formulated with t^, = to = and for the numer- 
ical computations t{ is chosen "large enough" for x{t{) to be close to the attractor, 
the chemical equilibrium point. The numerical value = is arbitrary as the ODE 

is autonomous. 

In contrast to this "forward formulation" , in the present work additionally the 
"backward formulation" = t» = is used. In fact, this is the more natural formu- 
lation for the identification of a trajectory on the slow invariant manifold which stays 
on this manifold during backward time evolution. However, the solution of the back- 
ward problem is much more challenging numerically since it is highly unstable and 
ill-conditioned for a dissipative dynamical system. We will refer to the first case with 
to = ti. = as forward mode and to the latter (if ~ t^ = 0) as reverse mode. Both 
modes can be seen as special cases of the general formulation (2.1). We deal with the 
numerical instability of the reverse mode by a collocation approach (see Section 4.1) 
with a fine discretization grid for the objective functional (2.1a) and the differential 
equation constraint (2.1b) and apply robust interior point optimization methods [31] 
to solve the resulting high-dimensional nonlinear programming problem (NLP). 



^x) := \\Jf{x) fh 



dx dx dx 
dt dx dt 



Jfix)-f. 




(2.2) 
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3. Theoretical Results. In this section theoretical results for the solution of 

the reverse mode problem formulation are presented. For a general two-dimensional 
linear system with distinct negative real eigenvalues and the nonlinear Davis-Skodje 
test model [6, 29], it is shown that for infinite time horizon of the variational problem 
the exact slow manifold is identified by the solution of the previously introduced 

variational problem. 

3.1. Linear Model. We consider the two-dimensional linear model 

2/iW = -A2/i(t) 
2/2(<) = (-A-7)j/2W 

with two time scales 0(A) and 0{X + 7) where 7 > measures the spectral gap 
(stiffness) of the system. In order to allow for a parameterization of the slow eigenspace 
by both state variables yi and y2, we apply an orthogonal transformation via rotation 
matrices R. 

R 

Hence system (3.1) is transformed to a; = Ax with 

^2 2 
T — X - 

2 ^2 



and the slow eigenspace is the first bisectrix x\ = X2- Since orthogonal transformations 

purely rotate the phase portrait of the dynamical system, the following considerations 
capture the general linear two-dimensional case. 

Theorem 3.1. Let x ~ Ax be a two-dimensional linear m,odel, A as in (3.2) with 
distinct (real-valued) eigenvalues —A and — (A + 7), 7 € M"*", fast and slow eigenspaces 
Af and Ag corresponding to — (A-I-7) and —A, respectively. Let x* be the optimal solu- 
tion of (2.1) witht^ =ti €R, g = 0, f{x) = Ax, and^{x{t)) = \\Jf{x{t)) f{x{t))\\l = 

\\AAx{t)\\l. 

Then for a// 7 > and tQ < tf it holds 

lim d{x*{ti),K) = lim inf \\x*{ti) - b\\2 = 0. 

to— ► — 00 to— ► — 00 beAs 

Proof. We assume w.l.o.g. the second variable being the progress variable, i.e. 
-ffix = {2}, and A = 1. The objective criterion $(a;(t)) can be computed as 

\\AAx{t)\\l ={xi{t)f (1 + 27 + 37^ + 27' + y 

+ (ar2(i))^(l + 27 + 37^ + 27^+^) (^"^^ 

+ xi(<) X2{t) (-47 - 672 - 47'' - 7*) . 

The general solution of the ODE x — Ax is 

:ri(t) = cie~* + C2e(-i-'^)* (3.4a) 
X2{t) = cie"* - C2e(-i-'')*. (3.4b) 
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Solution (3.4) is substituted into criterion (3.3) and integration over time yields the 
objective functional 



/ ' II AA:c(t) dt= r [2c?e-2* + (2 + 87 + 127^ + 87^ + 27^) c^e 

J to J to 



4^ ^2„(-l-7)2t 



dt 



= cf (e-2*o - e-2*f) - (e(-i- 



7)2fo _ e(-l-7)2tf j ^ (3 5) 



with ^ = 2+87+127^+87 +27 ^ expression Ci(c2) for ci as a function of C2 can 

be computed from (3.4b) which only dependents on C2 because of the fixed final value 
of X2{tf): 



X2iti) = cie-*f - C2e(-i-^)*f 



Cl(c2) = 



X2 {t{) + C2e(-i-T)*f 



This formula can be used to eliminate Ci from (3.5) leading to an expression h{c2) 
only depending on C2 (and to, U, 7, which are assumed to be fixed at the moment) 



h{c2) :-- 



(x;,') e-2*o e(-i-7)2t,g-2to ^ 24e(-i-'^)*fe-2*o 



-2tt 



+ 



-2it 



-2tt 



C2 



g(-i-7)2t, 2 _ 2x*'e(-i-T)*fc2 



_^e(-l-^)2to^2^^g(-l-7)2tf^2^ 



which should be minimal for identification of the optimal 02- The first order necessary 
condition for a minimum ^^^^ = gives a solution 



C2 



a;*fe(-i-T)*f - 4e(i-^)*fe-2*o 



e-27tte-2to _ ^e(-i-7)2to + (^ - l)e(-i-7)2tf ' 
Checking the second order sufficient conditions 



d4 



2e-27tt (e-2to _ e-2tf) + 3^ (e-2*fe-2^*^ - e-2*''e-2^*<') > Vc2, if > io 



guarantees 62 being a minimum. 

The solution £2 and Ci(c2) are substituted in (3.4a) evaluated at fixed final time 
tf yielding an expression for Xi{tf) additionally depending on 7 and to 



^i(tf) = ci(c2)e-*' + C2e(-i-^)*' = ^^±f?£_ 



1 + 



2e(-i-7)2tf _ 2e-27«fe-2to 
e-27«fe-2to _ ^e(-i-7)2to -|- (^ _ 1) e(-i-7)2tt 



2e(-i-7)2tf 

e-27«fe-2to _ ^e(-i-7)2to + (^ _ 1) e(-i-7)2tf 
2e-2T«t 

~e-27*t - ^6-27*0 + (^ _ 1) e(-l-7)2tfe2to 



(3.6) 
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with error term x quantifying the deviation from the slow eigenspace x\ = Xi. Finally 
in the limit ^ — oo it can be seen that 

lim X\(ti) = X2 

to — * — oo 

meaning the slow eigenspace Xi{t) = X2{t) is identified by a solution of the optimiza- 
tion problem. □ 
In Figure 3.1 the error term x is plotted. It illustrates that for increasing spectral 
gap 7 and increasing time intervals [to , t{] the approximation of the SIM improves 
while the error decreases exponentially. 




Fig. 3.1. The error term x (3-6) plotted against tg and 7 with t( = 0. 



3.2. Davis-Skodje Test Problem. The Davis-Skodje model (3.7) [6, 29] is 

widely used for analysis and performance tests of model reduction techniques supposed 
to identify slow invariant manifolds, 



dxi 
"dT 

dX2 

"dT 



= -xi 



-7x2 



(7 - l)xi + jxl 



(3.7a) 
(3.7b) 



where 7 > 1 is a measure for the spectral gap (stiffness) of the system. Typically 
model reduction algorithms show a good performance for large values of 7, which rep- 
resent a large time scale separation. Small values of 7 impose a significantly harder 
challenge on the computation of the slow invariant manifold. For reasons of adjustable 
time scale separation and analytically computable SIM, the Davis Skodje model is 
widely used for testing numerical model reduction approaches. We provide analytical 
and numerical results for the variational approach with the Davis-Skodje model. 
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Theorem 3.2. Let x = f{x) he the Davis-Skodje model (3.7), the slow invariant 
manifold defined by Ag := {{xi,X2) gM.^ \ X2 = o-'^d x* the optimal solution of 

(2.1) with t^=tf£ K, g = 0, /fix = {1}, and<^{x{t)) = ||J/ {x{t)) f {x{t))\\l. Then 
for all 7 > 1, a;^^ > 1 and to < tf holds 

lim d{x*{t{),As) = lim inf \\x*{t{) - b\\2 = 0. 

to— ► — oo to— > — oo 6gAs 

Proof. The Jacobian of / is given by 

/ -1 0\ 

Jf{x{t)) = (l+7)xi(t)+7-l _^ • 

V (i+xi(t))=* Ty 

$ {x{t)) in the objective function can be computed expUcitly as 

$ {x{t)) = {x,{t)f + 7^ {X2{t)f + ^"'^'^j' 



+ 



(-2-27^ + 47^) (a;i(t))"(l + 27^ + 67^) 
{l + Mt)f {l + xi{t)f 

{x,it)f{2^^+4^^) , 



{^ + x^{t)f {l + x,{t)f 



(2-f'x,{t) - 272 + 274 + 47^X1 (f) + 274 {x,{t)f) 

- Xl{t)X2{t)^ 3 ^. 

{1 + Xi{t)) 

An analytical solution of model (3.7) will be computed in the following. The first 
differential equation yields xi{t) = cie~* as a general solution. Equation (3.7b) is a 
inhomogeneous first order linear ordinary differential equation and the ansatz of the 
method of variation of parameters gives x-zit) = .^;2,hom(i) + ■'3^2,part(^)- The homoge- 
neous equations are solved by X2^hom{t) = C2C"'^', because X2,hom(^) — — 7C2C^'''* = 
—'yx2,honi{t). To determine a;2,part(i) the relation 

e-*C2(t) = (7-l)^^^-* + 7efe-^* 
(1 + cie-*) 

has to be solved for C2: 

, . / ^ , , (7 - 1) cie"* + 7c?e"^* , cie^* 
C2{t) = I C2it) dt = I ^ di = 



,{t)= J C2{t) 

Therefore, the missing part is 



e-T* (1 + cie-ty Cl + e* 



X2,pa.Tt{t) = e ^ = 



Cl + e* Cl + e* 
and the full solution of the ODE is given by 

xi{t) = cie"* (3.8a) 

ar2(i) = C2e-^*+^^. (3.8b) 
Cl + e^ 
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Now the criterion $ is integrated over time using (3.8). In the next formula r,, 
i = 1,2 represents a "rest" all terms independent of X2{t), hence also independent 
of C2, which annihilate after differentiation with respect to C2 afterwards. Making 
use of ci = xi{t{) e**^ = x^e'^'^ due to (3.8a) yields as an expression for the objective 
function only depending on C2 



h{c2) = n+ [\HMt)f dt 

J to 



= r-2 + c^(2'^V^'*°-2^'""''*0 



+ C2 

J to 



ci + e* ' ' (1 + cie-*)3 



dt 



-Mt) 

dh(c2) 



The necessary first-order condition for a minimum is applied. Setting — = 
results in an optimal 



C2 



-It>it)dt 

/ySg— 27*0 /y3g— 27tf ' 



The second order check 



= -y3g-27to _ ^3g-27tf q y^^^ to < tf, 7 > 1 

assures C2 being a minimum. 

An expression for X2{t[) can be derived by substituting ci and C2 in (3.8b): 

•^2 y^i) ^3e-27to _ ^3e-27tt ^ a;*!* + 1 ' 

The proof for the relation 

lim o ol \ „, e-T*^=0 (3.9) 

to— ►-00 -^32-27*0 _ -y3e~''')'*f 

will be given in the following Lemma 3.3. Because of (3.9) it holds 

lim X2{ti) 



*0 — oo Xi +1' 



which is the analytic expression for the slow invariant manifold of the Davis-Skodje 
system (see [6]). This completes the proof. □ 

Lemma 3.3. Under the conditions of Theorem 3.2 equation (3.9) holds. 
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Proof. In the following (x)^ is the Pochhammer symbol which is defined as := 
X {x + 1) ■ ■ ■ {x + n — 1). For the proof the integral in (3.9) is evaluated using Wolfram 
Mathematica® 7: 



/" 

J to 



(p{t) dt ■■ 



-27^e 



.2p(l-7)f 



(Ci7 + Ci + 76*) 



(ci+e*)' 



= :Vi(t) 



^ (2-7)„n! ^ ' *^ c?' 



n=0 



ci (7 - 1) 



=:V'2(t) 



With this result the values of tpi{t), i = 1,2 at t = to,ti are substituted replacing the 
integral in equation (3.9) and four summands can be regarded separately in the limit: 



*i,t,(io) 
*i,to(*o) 

*2,tf(io) 

*2,to(io) 



27^e 



2p(l-7)'t ('r'fo*f 



*fptt 



x^'c'-'7 + a;^'e 



7e*f) 1 



-7*f 



(a;*/e*t + e*t) (726-27*0 _ ^3g-27tf) 
2^2g(i-7)to (a;tfetf-y + a;*fe*f + 7e*°) e" 
(a;*/e*f + e*o)^ (-^33-27*0 _ ^3e-27tt) 



■7*f 



(X*/C*f (7 - 1)) (730-27*0 _ ^3c-27if ) 

(.?)"e-f 



(2-7)„"! ^ C.'fr."'f*^ 



n=0 



(x*/e*t (7 - 1)) (73e-27to _ ^3e-27tt) 



The limit of the first two expressions limt„^_oo ^'i,tf (io) = linito^-oc ^'i,to(^o) = 
is evident. The two terms 'i>2,tf{to) and \E'2,to(*o) contain a hypergeometric series. 
^2,tt(io) is absolutely convergent if x^-^ > 1, ^2,to(*o) is absolutely convergent if 
|e*''/x*''c*f I < 1 (generalized ratio test). This is always fullfilled for to small enough. 



Therefore, lim 



ta- 



*2,tf(in 



lim 



to- 



*2,t„(M =0. 



□ 



4. Numerical Results. We present numerical results for the examples investi- 
gated theoretically in the previous section. Additionally, numerical SIM computations 
for a simplified realistic hydrogen combustion mechanism axe shown. 

4.1. Numerical Methods. After suitable discretization the optimization prob- 
lem (2.1) can be solved as a standard nonlinear programming problem (NLP), for ex- 
ample via the sequential quadratic programming (SQP) method [25] or interior point 
(IP) methods, e.g. [8]. In particular, one has to decide how to treat the differential 
equation constraint and the objective functional. The easiest way is a decoupled it- 
erative approach, a full numerical integration of the ODE model with the current 
values of the variables subject to optimization. This procedure is called the sequen- 
tial (or single shooting) approach since it fully decouples simulation of the model and 
optimization. However, it is often beneficial to have an "all at once" approach that 
couples simulation and optimization via explicit discretization of the ODE constraint. 
This so-called simultaneous approach has the advantage of introducing more freedom 
into the optimization problem since the differential equation model does not have 
to be solved exactly in each iteration of the optimization algorithm. Especially for 
highly unstable ODE problems such as (2.1b) considered backwards in time a fully 
discrete collocation approach seems appropriate for the ODE constraint. On a prede- 
fined time grid the collocation method constructs polynomials obeying the difi^erential 
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equation at a certain number of nodes depending on its degree. For the numerical 
solutions presented in this work we use a Radau-method with hnear, quadratic, and 
cubic polynomials, respectively, [2]. 

The main difference between SQP and IP optimization methods for the solution 
of an NLP is the treatment of inequality constraints. Whereas SQP identifies the 
set of active constraints in the solution, IP formally couples the constraint violation 
to the objective function via a penalty term. Both methods finally use variants of 
Newton's method applied to the necessary optimality conditions, cf. [24]. For the 
numerical results presented in this work the NLP has been solved using the robust 
interior point method implemented in IPOPT [31] including linear algebra solvers of 
the HSL routines [13]. The required derivatives are computed using the open source 
automatic differentiation package CppAD [3]. Plots are generated using MATLAB®. 

4.2. Linear Model. Figure 4.1 and Figure 4.2 depict numerical solution results 
of problem (2.1) with the linear model (3.1) and small time scale separation 7 = 0.2 
and 7 1.0, respectively. Solutions for the forward mode and reverse mode are 
shown. In all cases X2 is chosen as reaction progress variable (parameterization of 
the SIM) and fixed at four different values: — 2.0, 1.5, 1.0, 0.5, for each of which 
the optimization problem is solved to obtain the coordinate of the second variable 
supposed to be located on the SIM (here slow eigenspace). The red curve is the SIM 
(slow eigenspace) which is given as the first bisectrix and the blue curves are the 
trajectories integrated numerically starting from those points (blue circles) that have 
been computed as solutions of the optimization problem. The red dot represents the 
equilibrium point (stable fixed point). Obviously the reverse mode gives solutions 
that are significantly closer to the SIM than the forward mode. 




X2 X2 
(a) Forward mode: = 0.0, tf = 10.0. (b) Reverse mode: to = — 21.0,tf = 0.0. 

Fig. 4.1. Results for the linear model (3.2) with 7 = 0.2, (a) forward mode: X2{to) = Xj", and 
(b) backward mode: X2(tf) = . 

4.3. Davis— Skodje Test Problem. Similar results are shown in Figure 4.3. 
Here the Davis-Skodje test problem is used for computations with forward mode 
(Fig. 4.3(a)) and reverse mode (Fig. 4.3(b)). In this case xi is chosen as reaction 
progress variable (SIM parameterization) and fixed at several values between 0.2 and 
2.0 for the computation of SIM points as solutions of the optimization problem. The 
spectral gap parameter is chosen as 7 = 1.2. 



A VARIATIONAL PRINCIPLE FOR COMPUTING SLOW MANIFOLDS 



11 





(a) Forward mode, to = 0.0, if = 10.0. (b) Reverse mode, to = — 8.0,£f = 0.0. 

Fig. 4.3. Results for the DavisSkodje test problem with (a) forward mode: xi (to) = a^*", 
and (b) reverse mode: x\ (t{) = Xj*. The red curve represents the analytically calculated SIM, the 
blue curves are the trajectories integrated numerically from those points that are the solutions of the 
optimization problem (blue circles), the red dot represent the chemical equilibrium point. 



Forward mode solutions show a larger deviation from the SIM (slow eigenspace) 
and a lack of invariance, whereas reverse mode solutions are highly accurate repre- 
sentations of the SIM. 

4.4. Simplified Realistic Mechanism. As a third example numerical results 
for a simplified realistic mechanism for hydrogen combustion are presented. The cor- 
responding full mechanism was originally published as a detailed hydrogen combustion 
mechanism by Li et al. in [20]. Ren et al. simplified the mechanism and used it for 
testing their ICE-PIC model reduction method in [28]. We use an adapted version of 
the simplified one. It consists of six chemical species (including the inert gas N2) and 
twelve chemical reactions as given in Table 4.1. Element mass conservation relations 
(in the general problem formulation equality constraints (2.1c)) for this mechanism 
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are 

xn + 2xB_2 + ^OH + 2.XH2O = 0.15 
a^OH + Xo + XHaO = 0.05 
2a;N2 = 1-6. 



Table 4.1 

Adapted version of the simplified mechanism of [28]. Rate eoefficients k are eomputed in de- 
pendance of temperature T as k = AT^cxp {—Ea/RT), where R is the universal gas constant. In the 
mechanism M represents a third body being any species with collision efficiency /h = 1, /h2 = 2.5, 
/oh = 1, /o = 1, /H2O = 12, and /n2 = 1- 



Reaction 

O + H2 
H + OH 
H2 + OH 

H2O + H 
O + H2O 
2 OH 

H2 + M 
2H + M 
+ H + M 
OH + M 
H + OH + M 
H2O + M 





A 1 cm, mol, s 


h 


P / kj 
-^a / mol 


H + OH 


5.08 X 10°"^ 


2.7 


26.3 


O + H2 


2.24 X 10°'' 


2.7 


18.5 


H2O + H 


2.16 X 10°^ 


1.5 


14.4 


H2 + OH 


9.62 X 10°^ 


1.5 


77.7 


2 OH 


2.97 X 10°6 


2.0 


56.1 


O + H2O 


2.94 X 10°5 


2.0 


-15.1 


2H + M 


4.58 X 10^^ 


-1.4 


436.7 


H2 +M 


1.18 X 10^9 


-1.4 


0.7 


OH + M 


4.71 X 10^8 


-1.0 


0.0 


+ H + M 


8.07 X 10^^ 


-1.0 


428.2 


H2O + M 


3.80 X 10^2 


-2.0 


0.0 


II + on + M 


(i.-lT X lO-'' 


-2.0 


499.4 



In Figure 4.4 results for the computation of a one-dimensional slow invariant 
manifold for the hydrogen combustion mechanism are shown. Solutions of the opti- 
mization problem (2.1) have been computed using the reverse m,ode. Again the red 
dot represents the chemical equilibrium and the progress variable xh^q has been fixed 
at different values between 0.0005 and 0.0180. The blue circles are the final values 
x{U) of the solution trajectories of the optimization problem and blue curves are the 
trajectories integrated numerically from those values forward to equilibrium. They ac- 
curately approximate the SIM; convergence of trajectories (dashed red curves) started 
from arbitrary initial values (red circles) to the computed SIM is visualized in Figure 
4.4. 

Figure 4.5 shows a two-dimensional manifold computed with the reverse mode. 
Two reaction progress variables xii^O and are fixed and the slow invariant mani- 
fold is approximated on a two-dimensional grid as a solution of a family of optimization 
problems. 

Acknowledgments. The authors thank Dr. Mario Mommer (IWR, Heidelberg) 
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Fig. 4.5. Results for two-dimensional SIM of the simplified combustion mechanism computed 
with reverse mode and aiHjoC^f) '^^'^ ^H2(*f) chosen as reaction progress variables, to = —5.0 X 
10~^,tf = 0.0, constant temperature T = 3000 K. The same arbitrary trajectories as in Fig. 4-4 ''^e 
shown in red. 
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